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ABSTRACT 


The finite element method is applied to the 
two-dimensional, oinvi Se wad compressible radial 
equilibrium equation foe axial CONDIeESSOTS. 
Isoparametric elements are used along with three-point 
GaussSlan integration for stiffness matrix evaluation. 
The ~radial equilibriua eqiation Vs pale 0 
quaSi-harmonic form for stream function formulation 
and results are presented using an isentropic flow 
assumption. Axial velocity profiles at rotor and 
Stator blade edges are compared with published 
performance data of the NASA Task-1 stage transonic 
compressor and with numerical finite element results 


of Hirsch and wWwarzee. 
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Ls INTRODUCTION 


A. PROBLEM STATEMENT AND OBJECTIVE 


The prediction of meridional flows within turbomachines, 
be they compressors or turpines, is a difficult but 
important part of the design process. The dirficulty arises 
from the presence of three- dimensional and viscous effects 
Within all turbomachines and the importance arises from the 


necessity to design accurately and efficiently. 


To simplify the problem of viscous, three-dimensional 
analysis, wWu [Ref.1] showed that this complicated flow may 
be analyzed by solving two interrelated flows: on? being the 
blade-to-blade flow describing the flow between rotating 
blades and the other being the meridional through flow which 
descrinpes the radial equilibrium. These flows are depicted 
in Fig 1. In addition, an inviscid and axi symmetric 
assumption is made in the through-flow thereby simplifying 
the flow to a two-dinensional, axi syametric, inviscid, and 


compressible analysis. 


Three methods may be found in current reports regarding 
the solution of the radial equilibrium equation. The first 
two are the streamline curvature method [Ref.2,3,and 4] and 
the matrix method [Ref.5 and 6] which 1s pasSically a finite 
difference technique. The third, a relatively new method, is 
the finite element method. As shown by Hirsch and Warzee 
fRef.7]°, the solution of the radial equilibrium equation by 


the finite element method is achieved by arranging the 


equation for the stream function in quasi-harmonics form. 


Due to the excellent results reported in Ref.7 and to 
further the research effort for finite eslement techniques in 
fluid flow problems, the purpose of this thesis is two fold. 
Pirstly, the goal was to formulate a computer program for 
solution of the radial equilibrium equation paralleling the 
steps as presented by Hirsch and Warzee. Secondly, after 
Suitable verification of computer results with those of 
Hirsch and Warzee, the goal was to compare computer 
predicted flows with measured performance data of the Naval 


Post Graduate School's transonic compressor. 


The purpose of this paper 1s t9 present a report on the 
results obtained thus far. In Section II, the derivation of 
the radial equilibrium equation is presented followed by the 
application of the finite element nethod to this equation. 
Section III describes the computer program in some detail. 
Section IV contains selected test cases which wer2 used in 
program testing and checking. In Section V, conclusions are 
presented along with recommendations for further study and 
WOrxK. Of the project: The appendices contain the prograa 
listing along with a sample test case for reference py the 
user. In addition, a list of references 1S contained for 


further reading on the subject of this paper. 
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Figure 1 - MERIDIONAL AND BLADE-TO-BLADE PLANES 





Lis “DaEORT 


Boe Hk DEREVATION OF THE RADIAL SQUILIBRIUM EQUATION 


The foliowing discussion is taken ‘from Ref. 7 with 
slight changes in notation. The basic turbomachin23 geometry 
to be analyzed is depicted in Fig 2. Although the machine 
noted is one stage of a compressor, a Similar analysis to 
the one that follows may be applied to other machines such 


ms axial turbines and mixed-flow machines. 


One begins with the Euler equation assuming the viscous 


forces to be negligible. 


sf ; F ¥) - : vel ete 


The continuity eguation, assuming unsteady flow is, 


at 4 v (pv) = O (II. A. 2) 


The First Law of thermodynamics in a fluid fiele 


yecomes, 


TYs= vh-VPlf (ez Tee) 
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AXIS OF SYMMETRY 


Figure 2 - TURBOMACHINE GEOMETRY 
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Substituting equation (Qe ae) into equation (EEA |) 


leads to the Crocco equation, 


—? 
W Vs (vx V) = Tiwi = vi (II. A. 4) 


where H is the total enthalpy. 


Assuming a steady and adiabatic flow, the energy 


equation becomes Simply, 


(¥-v)H= O Garg eS) 


which shows that along 2 streamline in a stationary syste, 


the total enthalpy is constant. 


In a celative syste, such oe the case in a rotor blade 
mow, tne total relative yeloOeGLey, q, can be expressed in the 


following form, 


—? a 
So yeoak = Ya (ieteenerc)) 


i _ 
where Ww is the constant angular yelocity and JU is the 


constant peripheral speed of the relative syste. 


Now, the Crocco equation in a relative system becomes, 


eee, 2 202 
sw - WW .(vxw} -Tys- Vh+ = = we (II.A-7) 


Parallel to equation (EEeh. 5). foe the stationary systea, the 
energy equation, assuming steady and adiabatic (relative) 


flow in a relative systen, becomes 


4 


12 





_ (II. A.8) 
(Ww v ) He =O 


where H_is the relative total enthalpy expresed as follows, 





Ho =h+ Ww uy? (TTA 9) 
- “9 


From the following velocity diagran, 


13 


Be 






equation (II.A.9) may be arranged as follows. 


Since, 


3 Z 2 2 = (PIA 10) 
Win t Wg - W = +(Vop =) 


then, 


2 


We Vorn Vo ~ 2UVo a (Tar A. TY) 


and, 


2 


Wie ve le = ls (EA 2.2) 


SUDStLitUfIG 2cliation (LieA.l2), imtoo €@quaeton (2724.59) 


leads to the following relation, 


He=h + ¥ - UNp = H- UVe Caron 13) 


BQUat1on {2E2A. 8) Shows that He is constant aloag a 


streamline in a reiative system. 


Upon integrating equation (11.A.8) between the rotor 
inlet and outlet, the Buler equation for turbonachines is 


Found, 


Our 


il 


An | A (iu Ys)} 


[SJ iN 


(II.A. 14) 
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It may be shown [Reft.9] that by circumferentially 
averaging equation (II.A.1), and under the axi symmetric 


flow assumption the following relation is valid, 
=? = . 7 
-V«(V4V \s T-¥S - VH + Fy te) (IT.A.15) 


where PL is the body force of the blades acting on the fluid 
and all variables are mean values along the direction of the 
Circumference. Hence equation (IT.A.15) iS an approximation 
for axi symmetric flow. As a final note on equation 
(II.A.15), Since the viscous forces were neglected in 
equation (II.A.1), there must be a force introducing the 
entropy variations along the blade. big aelrs force is 
proportional to the pressure loss coefficient and is labeled 
Ep, Che dissipative force. F&F, produces work Which an turn 
produces entropy production radially along the blade. Under 
the axi symmetric assumption, entropy varies axially and 
radially only and is assumed to be proportional to the 


pressure loss coefficients {[Ref. 2 and 8}. 


- Due to boundary conditions imposed on the problem and 
the axl symmetric assumption, Sylinde Leal 
coordinates, (c,@,Z),will be used in all subseguent analysis. 
Phererte es, equation -(1E.4.15), 20  ~cylandrical coordinates 


and with axial symmetry is as follows, 


|} 


Ve 2 ). HM = 
2 , (Eve) Vl $5Ve - Vs) - oe Ce He ie ——- 


q- 


3S (Eve) + _ =, (eve) = Fe ee) 





Oo Oo 


. > sf — 3S 
Ve ($a \e - 4) - Yes (ave) = <u { 32 -F, (IT.A.18) 


It 1S important to note here that under the axi 
Symmetric assumptions, eguation (II.A.15) reduces to the 


following, 


— yp ad 
V- fos O (Tren. 19) 


Likewise ina relative system (rotor), the axl symmetric 


assumption leads to the following, 


Equation (II.A.16) describes the neridional through flow 
radial equilibrium equation for the finit2 element method. 
Since one 1S concerned with the meridional planes, the 


following derivative expression is taken from Fig 3. 


a a oi: (ITI A.21) 
Vina Ve or + VE oS 


Therefore equation (II.A.17) reduces to, 
d 
RF 5 7 Vim Syn (KR Vg) (Tie Aa2) 


which reveals that ina duct where there are no blades and 


therefore no blade forces, angular momentum is constant 


16 





aa a OO OO 


Be 


along a streamline. In that case, 


a 
5m (RVo) =O (II .A.23) 
As Shown in Ref.9, the circumferentially averaged 


continuity equation is the following, 


52 (xe bVe) sf x, (pb Vs) =O (II.a.24) 


where bp is the blockage factor defined by Hirsch and Warzee 
as the tangential area reduction due to the thickness of the 


blade. 


c 
b= i-s (Prone 5) 


where t is blade thickness and s is blade spacing. 


VW 











MERIDIONAL PLANE 


Figure 3 - 
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One further step in the formulation of the radial 
equilibrium equation for solution by the finite element 
method involves introducing the strean EUHe LION. In 
cylindrical coordinates, the stream functions are defined as 


follows, 


jo 3 
Vo = kb oe 


(ihe 6) 


oY (teen) 
ot 


Substituting these expressions into equation (II.A. 16), 


the equation becomes, 
bola Bi 1 (at - 
= Vou 
-B 5e [RVe)- FF 


The right hand side of equation (II.A.28) is applicable 


(TI.A.28) 


to the absolute flows in the stator and duct regions. FOr 
relative flows such as those in the rotor, the right hand 
Side is modified by replacing the total enthalpy,H, by the 


relative total enthalpy,H and the quantity Ve/R is 


Q v 


replaced by We/R. 


AS a last assumption in the formulation of the governing 
relation for the meridional through flow radial 2quilibrium 
eguation, both the radial component of the body force, FL 
and the radial component of the dissipative force,F,, are 
hegiected. This assumption, [Ref.1,8] does not hamper the 
accuracy of the results for conditions at design speed. 


Even though oublished compressor performance data used for 


19 


the test case in this thesis was obtained at 9.5 design 
speed, these force terms were also neglected in the computer 
program. As will be shown later, this assumption could 
possibly have had adverse effects on the predicted axial 


velocity profiles at the rotor hub and tip regions. 


The final representation of the meridional radial 
equilibrium equation to be solved by the finite element 


method 1s as follows, 


ad [dv a {,,a¥ 
(KS eal 3a (ke) + f - oO CT Teste 9) 


where, 
K = ~&b (II.A.30) 


and 


As 26 


Be VUGs FINTTESEbetety MET OD “APPLIED FO THe RAPEAL 
EQUILIBRIUM EQUATION 


In order to formulate equations (II.A.29) through 
(II.A.31) in matrix form for solution by the finite element 
method, one must apply a weighted residual technique to the 
equations for numerical solution. The weighted residual 
method used here is the Galerkin's Method. The following 
discussion is taken from Ref. 7 with only slight changes in 


potation. 


REWr1ting. G@quation (Li. A.29)" ahd divicing througa by k, 


one has, 


20 





OB Be re OO 


Tn 


(its Be 1) 
ae (ee) 42 = 
eg poe 3e/t 33 

where this equation represents the flow in the volume,V. 


The boundary condition for this partial differential 


equation, after dividing through by R, is, 


zk — (4, ) =O (II. B. 2) 


Where this egquatioa solves the flow on the closed boundary 


6, the vVoluge,or,s. 


By applying the weighted residual process to equations 
(Li. 8.1) and -{El.3,.2) (and Using an arbitrary weilgatinag 


FUNCTION, Wityz), one bas 
{ 
wee 2), IV e | wie) Cue dS =O (II. B. 3) 
V s 


where Ee and E coy are the volume and surface resijuvals 


respectively, or, 


Lyd ae d (34 
rae £ [Releet)« SGA) Leo ones 


| 
eon Eg ra + XM, (y) - , (II.B.5) 


he) the  SOLUcion. “tO eam@ation (io. 8. 1).was exact, pots Cot 


and r... would be equal to Zero. 
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(II.B.2), one may analyze the egyuation as follows. 


2a 


On the surface,S,where W 1s specified, 


W - y (II.B.6) 


and, 


HY, -7% oO (II. B.7) 


ov 


Similarly, on the surface, where n° °,Ss where 


ait 


( Faby, B. 3) 
q =O 


Sas, =O. US, Z e (£i./5s:9) 


Due to the axi symmetric assumption, the final equation 
Will not involve dV and dS but the intersection of dV and ds 
With the meridional plane. Therefore, one must transform the 
volume integral,dV, to a surface integral and the surface 


Integral,ds; €O 42 laine integral: 


Hence, Let, 


dQ = intersection of dV and meridional plane 

dC = intersection of 4S and meridional plane 
= dV= 27 Rd 

dS= 2@wedc 


With this transformation, one may rewrit2 equation 


ZZ 








(II.B.3) as follows, 
sv\ a/ ab) 
\-wieial alk | * 3 (def 20 ee 
{ w(2k Sa 
Cc 


where on the contour, ee. 


(II.B. 10) 


One must now integrate the first term in equation 
(1228.10) by parts to Obtain Gae following, 


-\w {Se (eS n+ (x) onda ~ Wwode 
fe eae bo. 7 wide andl =O 


Inspecting the first term in equation (I1.B.11), one may 


(or. 8. 11) 


use the following integral theoren to simplify further 


[aabd.n = [onde (ft -3 510) 
Gq 


SL 


Rewriting a (TISBsll) ,»dives, 
{Py YW. 3Y \W 
Swi xe "e+ Yn, dC- fda ‘ 2 se M) Je 


-\we 4 2ndl =O 


C 
Finally, since 


(it ori) 


yp a oY 
aA = se UR 7 "3 (i 


SGuation. (see 14) reduces to tie following, 
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ey. 4 


ait 


ste 


my). {Wlda=o {TI .B. 15) 


One now has the final equation in the form for use by 
the weighted residual method using any arbitrary weighting 
£UnCtiOnN, W(r;,zZ); As noted previously, the Galerkin's 
Method will be used here which implies that the weighting 
functions are the same functions used in approximating the 


stream funeeion, ¥ 


Before applying the finite element method, one must 
discretize the continuum and then approximate the unknown 
function iY, by a set of polynomials. For this particular 
problem, eight-noded iso-parametric elements were chosen for 
discretization, see Fig 4, and the following approximating 


functions were used. 


ps 
HM 


g 
2 
l 


= 


Ne (G0) % (II.3. 16) 


where, 


shape functions 


Ni \ 1) 
¥, 
u 


location within the element. The shape functions, 4: , 


here are defined by the following relations as shown in 
Refs 10: 


value of V at the node 


value of VY at med ae Dera 


used 
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Ni (8) = El 86:)CI+ ype) (EE: HM; - 1) 
Ni (§,7) = z (1- ) (+49) (i220) 
Ni (5) = (I+ 64: ) (| 4) 


where the following coordinate transformations are used, 


t= > N (4,4) Y; 


(II .B.18) 


5 
z= z Ni (4) 
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Figure 4 - ISOPARAMETRIC QUADRILATERAL ELEMENT 





At this point, one is ready to apply tke Galerkin Method 
to equation (II.B.15) by substituting equation (II.8.16) for 
the unknown function Y , and N; for the weight function, W, 


which yields, 
. g ; 

(else Salas). SD) da [honk da = Octo 
Pet 2 


This integration yields the folloving system of 


equations which is solved for the unknown nodal, 


Ei Rigs kin Y 5 





; (II.B.20) 
Kui Kww Yu s, 
where, 
ANC INS ONS IN: 
Kis I, iS pe ~ ye a d ee ete eB 
St 
and, 


I. = VEN d2 (Pi se 225 


In addition since both the ‘stiffness matrix',&, and the 
right hand side vector, [FJ], are functions of Y , the systen 
as defined by eqrations (I1I.8.20) through (II.8.21) must be 


solved iteratively. 


At this point one has the total finite element 
formulation of the radial equilibrium equation as defined by 
SGuations (fl, 8.19) “end (iicowc0). The problems which 


remain to be clarified ar@ basically two fold. Firstly one 


2/7 





must evaluate the integrals in equations (II.5.19) and 
(IZ.B.2) by numerical methods, and secondly, the solution 
procedure for the non-linearity must be formulated. In Part 


C, both of these final steps are presented. 


C.  NUMERTCALSINTEGRATION OF STIFFNESS MATRELC AND SOLUTION 
PROCEDURE 


1. Numerical integration of the stiffness matrix 


As noted in Section II.B, 2valuation of equation 
(II.B.21) must be perfocmed numerically. In addition, one 
realizes that the derivative exoressions enclosed within the 
interval must be evaluated by a coordinate transformation. 


This is done in the following way, 


Since, 
3 
Y= 2 Ny (s 1) t 
(Tene. Is 
Ea 
z= 2 Ni (Sy) 2 
then, 
SN; | ANG Je | SONU Or 
ye dE 38 yr} 
" 7 (Tire C..2) 








and 2 Matrix fora, 


id) [Sear | (an 
. 3 Wel) 

4 = " i (II.C. 3) 
ONG gt ry eae 


my VE Aer 


Furthermore, defining the Jacobian matrix as, 


3¢ dr 
a 
J = 7 " (Tis. 6) 
Ze Sr 
1% 
then 9 DY  @ividing both Sides of equet ton {11.C.3) by J, one 
has the following transformation, 
AN fam 
st zy +6 
ee rr (ks C..5) 
; (s| 
oN: ONG 
av | | 3% 


In addition, it has been shown [Ref.9] that 


Atdr | s| dg dy (II.C.6) 


NOW, With. egcuations § (11.0.5) and (Ue Cay equation 


(II.B.21) becomes the following, 


Zo 









Ki > (, se ey Yor ¢ Ly athe (TT eed) 


Equation (Lie ss 72) is best an using the 
Gauss~Legendre integration method since it is of the 


following f£o0rn, 


ee = |e md§dy (21ers) 


Or finally, [| kerst0}],; 


z 148. 5 (8, We (Ttnc.9) 


where A.and B.are Goeifieztents (Fig 5) £or both “wo and 
b 


three point GausSian Quadrature. 


At this point, one has the tools to calculate all 
the elements of the stifiness matrix. In like manner, the 
right hand side vector,f, 1s calcualted by numerical 


Piteqrataon. 
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NUMBER OF 


GAUSS IAN qi 


POINTS 


e577 2207 F 1.00000 00000 


0.77459 66692 O45 59 5035555 


0.00000 O0000 0.88888 88888 





Figure 5 - GAUSSIAN INTEGRATION POINTS 
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2. Solution procedure 


The following is a synopsis of the basic solution 
process. Specific details concerning equations and methods 
of computer coding are covered in the proceeding section. 
The proceeding 1S meant t9 give the reader a preview of the 


solution process. 
an DESsecretization 


Initially the machine under analysis is 
discretized into eight-node iso-parametric elements. The 
axial calculation stations are placed arbitrarily in the 
duct regions and along blade edges and centers for the rotor 
and stator as shown in Fig 6. At this point the syst2n 


topology and nodal coordinates are specified. 
Db. Initialization 


To begin the iteration process, one aust assume 
an initial internal stream function, velocity, and density 
distribution. In the program, the initial internal stream 
function was assumed to be that of the outer boundary 
throughout while the velocity and density distribution was 


Vv 
assumed to be that of the inlet. 


a2 
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c. Calculation of thermodynamic variables 


Before calculating the right-hand side vector,f 
, one must obtain distributions of angular monentun, 
enthalpy, and entropy. This is done by first calculating 
the thermodynamic variables at tne inlet axial station from 
the given inlet conditions. In order to proceed axially 
through the machine to calculate the nodal angular momentun, 
enthalpy, and entropy, tne following three equations derived 


in Section III are used. 


H = Cot = constant along a stator streamline 
Ho= C a ae) = constant along a rotor streamline 
R- pls z g 
Vo = constant along a duct streamline 


An example of this calculation procedure for the duct region 


is shown graphically in Pig 7. 
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Figure 7 - DUCT ELEMENT 
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In this figure, the angular momentum at point X 


is equal to the angular momentum at point Y. More formally, 


5 
Vale ENG), = WVp)y SEER 


Since the previous axial station's thermodynamic 
“variables are known, one must now find the values of g and 
yn at point Y. This is done iteratively in the following 


way. Since 
5 
Y= P| = 2 Ni (@.4) ¥: (Tie 4222) 
U=3 
and along the left side of the slement, 
_ -4 (II.C. 2. 3) 


then equation (II.C.2.2) aav be solved for | by a Suitable 
lteration method. As will be shown in the next section, a 
halt-interval method was used to obtain the unknown y 

Once H 1S known, then equation (II.C.2.1) 1S solved for 
the angular momentum at point X. The rotor and stator are 
handled in a similar fashion. In addition, the rotor and 
Stator deviate the flow creating a three-dimensional flow 
field between the blades in the respective blade row. Low 
Speed cascade correlation data [Ref.13] was used to 
Calculate the effective turning angles in the rotor and 
stator. These effects are calculated beforehand with known 
mass flow rate and uniforn axial velocity assumptions at tne 
rotor inlet. The results of these calculations are part of 
the input data routine in the form of relative and absolutes 


flow angles at the rotor nodes and absolute flow angles at 


36 








the stator nodes. This will be shown more exactly in the 


next section. 


d. Calculate matrices 


At this point, the right hand side vector,f, and 


the stiffness matrix ,k, are calculated. 


e. Solve system of equations 


The system of equations as shown in equation 


(II.B.20) is solved for the nodal stream function. 


Eg Perrorml tTelaxiation iterat ron 


Due to the strong non-linear properties of the 
system of of eguations, the following iterative scheme is 


necessary. 


YV+\ y\ “~~ yi 
c c L C 
where @ is the under relaxation factor. As will be shown in 


Section III, this scheme is performed only in certain 
regions of the machine and in addition after a specified 


number of iterations. 


g. Update velocity and density profiles 


Using the current nodal distribution of the 
Stream function, axial and radial nodal velocity components 


are calculated along with a new nodal density distribution. 
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Again, this calculation procedure will be shown in the next 


section. 
h. Test for convergence of Y 


Stream function convergence criteria is now 
tested and Will determine if further iterations are 
necessary. The solution is said to converge if the following 


equation holds for alli nodes. 


Y\ 
a ai Y. 
oo 


t 


n+ | 


mie (iis 235) 


where € ius a designated requirement for convergence. 


i. Summary 


In summary, the eight steps involved in the 


solution are noted below; 
(1) Discretize the continuun. 


C2) Assume an initial stream function, velocity, 


and density solution. 


(3) Calculate the nodal thermodynamic variables 


from the given inlet conditions. 


(4)Form the right hand side vector, f(r,2), and 


the stiffness matrix, kK. 


(5) Solve the system of equations, given by, [kK] 


= {(F] for a new stream function distribution. 


36 






(6) Perform relaxation iteration if required. 


(7) Calculate new nodal velocity and density 


distributions from the current stream ftnection solution. 


(8) Test the solution for convergence, and if 
required, repeat steps (3) through (8) using the current 


nodal stream function values. 
This concludes the solution description and now 


one is ready to more completely understand the computer 


program which assembles the preceding eight steps. 
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Eri. THE PROGRAM 


A. OVERALL FLOWCHART AND DESCRIPTION 


The overall flowchart of tne program is depicted in Fig 
S, Those blocks denoted by the letter 'S' are subroutines, 
While the remaining calculations are an integral part of the 


Main program. 


After proper dimensioning of all arrays and subsequent 
initialization, the inout data are read and then printed. 
This not only presents a physical picture of the problem but 
also serves as across check to the user for correct data 
insertion. In addition, a subroutine is availabl2 to obtain 
a computer drawn plot of the mesh (Fig 6) and 1s a further 


check on proper data input. 


At this potnt ail the necessary variables nave been 
stored and the iteration counter for strean func Lon 
convergence is set. @#ith the current nodal values of Y and 
the given inlet thermodynamic conditions, the thermodynanic 
variables throughout the machine are calculated. From the 
calculated values of enthalpy and anguiar momentum, 
(isentropic flow is assuned), the right-hand side vector is 
Calculated followed by the stiffness matrix calculation 


(equation 17.8.21) . 
The system of equations (equation (I1I.B.20)) 1s row 


solved for the few node! Stream Function Gastriovtion. rr 


iS here Where fOr all 2terarions “Out “Ehe  Firse chat a4 
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relaxation factor is applied as noted previously in equation 
(Ese a2e4) The reasoning behind not applying the 
relaxation scheme to the value of nodal Y after the first 
1teration is the fact that the first iteration produced a 
close approximation to the correct strean function 
distribution. With this clos2 approximation to the strean 
function came a velocity and density distribution which in 
turn was near the correct solution. It was found that if 
the first iteration was relaxed, the second iteration became 
unstable since in fact the velocities and densities were 
themselves farther from the true values than were assumed 


Ligeuredt y. 


After testing the nodal stream function for convergence 
DY USC Of GGYation (Lllec.2.5), tne calculation precess is 
either repeated or ceased by virtue of convergence or 


limiting the number of iterations. 


As stated previously, low speed cascade correlation data 
fRef.13] were used to calculate turning angles in the blade 
regions. These angles were assumed constant throughout the 
Solution and not refined after subsequent Lterations. 
Further work on the somputer program could entail an 
additional computational routine which would calculate the 
new tient ad angles after eacn iteration. A sample 


calculation of rotor turning angles is shown in Appendix D. 


In the following sections the orogram structure is 


examined in more detail. 
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DIMENSION ARRAYS 










INITIALIZE ARRAYS 
READ IN DATA 
SET ITERATION COUNTER 


CALCULATE THERMODYNAMIC 





S 
VARIABLES FROM GIVEN 
CALCULATE RIGHT HAND 
SPEDE VECTOR FROM GIVEN S 






THERMODYNAMIC VARTABLES 





CALCULATE STIFFNESS MATRIX 
SOLVE SYSTEM OF EQUATIONS FOR y 


PERFORM RELAXATION 
ITERATION IF NOT Ist ITERATION 




















CALCULATE NEW VELOCITY AND DENSITY 
DISTRIBUTION FROM NEW STREAM FUNCTION 










STREAM FUNCTION 
CONVERGENCE SATISFIED OR 


ITERATION LIMIT 
REACHED 








YES 


PRINT FINGTS 
ELEMENT RESULTS 


Figure 8 - PROGRAM FLOWCHART 





INCREMENT ITERATION COUNTER 
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B. THE MAIN PROGRAM 


1. The input routine 


The following is a description of the input data 


required by the program. The data are arranged into 


categories described in the following manner. 


a. category 1 


Problem identification. 


ox category 2 


Number of nodes and number of elements. 


GC. “Cavegory 3 


Node numbers, nodal coordinates and 


blockage factor. 


d. category 4 


System topology. 


@®. category 5 
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twelve 


nodai 











function. 


Element type, duct, ©Otor, or Stator, 


category 6 


Absolute flow angles for rotor and stator nodes. 


category 7 


Relative flow angles for rotor nodes. 


category 8 


Inlet thermodynamic guantities. 


category 9 


Physical constants for fluid under observation. 


category 10 


PLTst esto tate Of Internal Stream  fuinpet 2on. 


category 11 


Node numbers and Sooe. . ec nodal stream 


category 12 
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Node numbers where the right hand side, f(r,Z), 


mS tO De Calcutated. 
Before describing in detail the format to be 
followed for data insertion, it is important to note the 


following assumptions. 


(1) Uniform flow conditions at inlet and 


Oowelet. 

(2) Uniform flow conditions at rotor inlet 
Lor Ca levitation of appropriate turning angles. This 
assumption is necessary to calculate the values of rotor and 


Stator flow angles. 


With “this “lf - @ina@ , ene discussion Will 


continue. 


The following describes each category in more 
d2tail. 


Category 1: 


Format: (20A4) 


Niiiper Of Caras: 1 


Procedure: Enter the title of the problem in 


columns 1-20. 


Category 2: 


ROEdare 21 10) 


Number of caris: Equal to the number of 


nodes in the system. 
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Procedure: Enter the number of nodes in 
columns 1-10, and the number of elements in 11-20. Both 
integers must be right justified. 

Category 3: 


Fotmate (110g 3F 10.0) 


Number of cards: Equal to the number of 


nodes in the system. 


Procedure: Each )Card Contains. tne nod 


M 


number followed by the Z coordinate, R coordinate, and nodal 
Dlockage factor. The coordinates are in dimensions of 
inches. 
Category 4: 
PoOrmrac: (91S) 


Number of cards: to the number of 


ty 
+Q 
ft 
fw 
-~ 


elements in the system. 


Procedure: sach card contains nine integers 
Pegne Justified in columns 5, 10, 15, 4tc., through 45. the 
first integer is the elenent number followed by the eight 
nodes associated with that element. It 1S important to note 
that the nodes ar2 read in starting with the upper fright 
hand node and proceeding in a counterclockwise fashion 


around the element. 


Category 5: 


Format: (2110) 


Number of cards: Sqgual to the number of 








elements. 


Procedure: Enter the element number in 
columns 1-10, followed by the integer '1! (ater) ta? 
{cotor), or '3' (stator) describing the element as either in 


SoMIdeG, .OLOL, -Or Stator rod Lon. 
Category 6: 
FOEMats (O%740,1 10,7 10.0) 


Number of cards: Equal to the number of 


rotor and stator nodes olus one 'STOP' card. 


Procedure:  2nter th] node Humper {right 
justified) in columns 11-20 followed by the value of the 
associated absolute flow angle in radians in columns 21-30. 
The last card in this category is a 'STOP' card entered in 


columns 7-10. 
Category 7: 
POEMeGt. Oat yO ye oe) 


Number of cards: Equal to the number of 


rotor nodes pilus one "Sree" card. 
Procedure: Enter the node number (right 
justified) in columns 11-20 followed by the value or the 


associated relative flow angle in radians in columns 21-39. 


THe last card in this Gategory isa 'S 70?" card. 
Category &: 


Porn ac: (7210.0) 7, (2 10.0) 
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Number of cards: 2 


Procedure: Enter the following quantities in 


the prescribed order and with the noted dimensions. 
First card 

Mass flow rate: (lbm/sec) 

Inlet axial velocity: (ft/sec) 
Outlet axial velocity: (ft/sec) 
Pnjlet total density: (lbmft’ ) 
Enlet Static density: onyee 5 
Inlet total ovressure: (1647 in ) 
Inlet total temperature: (°R) 


Second card 


Gateqony 3: 
Format: (3F10.0} 
NUMDer Of cards: 1 


Procedure: Enter the following quantities in 


the prescribed order. 


Gao Conctaut.. (12-lor7 ou &) 
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FO rrr 


Ratio of specific heats 


Constant pressure specific heat: 
(BTU/Lbm—-°R) 


Category 10: 
Formac: (.710'.4)) 
Number of cards: 1 


Procedur?: Enter the Tirst estimate of the 


1mteornal “Stream. LUNCtITOnN tO be WSSd im the EierSst 2 teoration. 
Category 11: 
FOrMaty (ox PAG, ? WO7F10.0) 


Number of cards: Equal to the number of 
nodes having a svecified value of the stream function plus a 


‘Ss LlOr* ‘Cara. 


Procedure: This set or cards allows the 
Stream function boundary conditions to be read in. A 
tyelcdl Card CONntalnS an Lateger, right Justifved tn colunns 
11-20 , which is the node number, followed by the vaiue of 
the specified stream function in columns 21-30. The last 


Gata 22 2. S10 card. 
CavedoL you. 
Pormat: (5X,A4,110) 


Number of cards: Equal to the number of 


nodes where the right hand side is to be specified. 
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Procedure: Enter the node number, right 
justified in columns 11-20, where the right hand side is to 
be calculated. Again, the last card in this category is a 
'sTOP' card. 


After all the data has been read by the program, 
the input data is printed and the mesh is plotted for 
verification by the user. The sample format is shown in 


Appendix C. 
This concludes the input routine. The next 


section describes the calculation of the stiffness matrix, 


Ke 


2. stiffness 


Ww) 
\3 
\f 
(ct 
in 
| 
>< 

(D 

< 

fp 
+- 

CG 

py 

a 

ts 

O 

an | 


as shown previously 1n Section ILE.C.1, the following 
equation describes each term in the eight by eight elemental 


act. t 


wie Se a \(ay ‘ let (Tldedy (eet 3). 24) 


In addition, 'k' is defined in the following way 1m OLrder. to 
numerically integrate th=> equation. 


1 


= as 


= PNGIG, i. 2 Ni (¢ y) 4b 


(77. Bi ese) 


where b is defined as the elemental blockage factor taken as 
an average over the eight nodes of the particular element 
and (FM) are the defined Zauss-Quadrature integration 


points. 
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=i) 















iS ae 


| FIND SHAPE FUNCTIONS AT 
GAUSS-QUADRATURE POINT ''I" 


FIND JACOBIAN AT "I", [J] 







FIND DENSITY AND "R" COORDINATE AT "I" 






FIND DETERMINANT OF JACOBIAN AT "T' 








N t t 
FINDS N, oN, Fam 


¥§ | 


FIND 











a tt | 









FIND ie Ano? 






FIND fra- 1 Tar "re 





PERFORM INTEGRATION AT ''I"! 


YES YES 









II G.T. NUMBER 
QE ELEMENTS 





NO 


ASSEMBLE STIFFNESS MATRIX 
W/O REGARD FOR BOUNDARY 
CONDITIONS 


STIFFNESS 
ee Ne MATRIX COMPLETE 


Figure 9 - STIFFNESS MATRIX EVALUATION 
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Fig 9 depicts the flowchart for both elemental 
stiffness matrix evaluation and the assemblage into the 
system stiffness matrix. More specifically, the figure 
ShowS a three-point Gaussian Quadrature scheme but can be 
Changed to a two-point scheme by simply integrating four 


times instead of nine as shown. 


The actual coding of the stiffness matrix evaluation 
and assemblage may be found in lines S$TRO3510 through 


STROU770 in the computer progran. 


3. Solution of systems of equations 


> ke ae ao cee —_— =~ —_—e ee ae ee = —e ee ee eee 


At this point, the system of eguations are modified 
for the boundary conditions and solved for the nodal strean 
function values. An equation solving routine, DS ENO; 
available in the system library was used for this purpose. 
It was found that no conparable savings was realised by 


using a banded equation solver. 


i Iteration schemes 


Ke hoted » previously in Section [1.C, 2 reeaxation 


scheme is necessary for convergence to a solution. 


Two ave tine: dirferences with regard to the 
iteration method were noted from that of Ref.7. Firstly, it 
was found that relaxation was necessary only in the rotor 
and stator elements and also in the duct region between the 
FOUOL OUtlek and stator aniet. Secondly, due to the extrene 
Mon linearity in the \rotor-stator areas, - 4° S¥itch Was 
required which changed the sign of 4 in equation (II.C. 2.4) 


as required for Stability of conVergerce. Ciarirticatton oF 
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this change follows: It was found that during the initial 
three or four iterations, the stream function values of the 
rotor-stator nodes sometimes exceeded the value of the upper 
boundary. Due to an absc2nce of sources Within the domain 
of solution, this occurrence was incompatible with the 
boundary conditions. At this point, it was necessary to 
make & negative in equation (II.C.2.4). During subsequent 
iterations, as the solution conversed, the rotor-stator 
regions became stable and the sign of & was returned to its 
positive value. This iteration proved to stapdilize the 
solution with respect to stream function values and 


Velocities, 


The iteration procedure is coded in tre comouter 
program from lines STROS5O070 through STR05200. 


Once convergence is obtained or the number of 
iterations have reached the limit imposed py the user, tne 
results are displayed. A sample output is shown in the 
Appendix. In addition, the units of all dependent variables 


are the same as those noted in the input routine. 


Ce THE SUBROULIUeS 


The following describes each of the six subroutines in 
the computer program. Each subsection contains a list of 
Calling arguements and for subroutines FPCAL, SLINE, and VEL, 
a basic flowchart. In addition, for those subroutines whose 
mathematical theory was not presented in Section Iff, a 


brief treatment is also given. 


=, 








1. Subroutine shape 


= ap a a 


This subroutine calculates the shape fuUnCcLions 
(equation (II.B.17)) at the values of & and i requested 
in the arguement list below. 


SUBROUTING SHAPE (5,272,527) 


E = value of § (input) 


N 
il 


value of y (input) 


SF = eight by one vector of the 2ight shape functions. 
2. Subroutine jacob 


JACOB calculates the Jacobian matrix as defined in 
equation(II.c.1.4) for the value of Bi denoted in the 


arguement list. 


SUBROUTINE JACOB (21,21,D,2,2C3$,ZC$,RJAC) 


E1 value ory (input) 
Z1 = value of 5 (input) 
Ne 
D) = Pight by one Vector of 3g (calculated) 


SNe 
E = eight by one vector of Sy (calculated) 


RC$ = eight by one vector of the 'r' coordinates of the 


nodes associated with the element (input) 


Z€5 =-eignt by One vector or tie 2" cocrdimates of (the 


nodes associated with the element (input) 


ROAC = two DY LeO Jaca Dilan Watrix (Output) 
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In addition, the subroutine assumes that the vectors 
RC$ and ZC$ contain element coordinates arranged in a 
counter clockwise fashion beginning with the upper right 


corner node. 


3. Subroutine sline 


This subroutine calculates the thermodynamic 
variables throughout the machine given the inlet conditions 
as described in Section II.C.2. The calling arguements are 


defined below. 
SUSRCUTLUE Shine (UCN EET RC pe olathe tel ole) leg 
NODE NNODET Celi RR, ALP NG, UNWELL 8S, HS) 
UINGET = Inlet axial velocity 
RC = Nodal 'r' coordinates vector 


Nodal stream function vector 


Pot 


iI 


WRL Nodal angular nomentum vector 


H = Nodal total enthalpy vector 


UVEL = Nodal axial velocity vector 

VVEL = Nodal radial velocity vector 

TVEL = Nodal absolute tangential velocity vector 

NODE = Matrix containing nodes associated with tne 
element 

INLET = Vector containing node numbers at inlet station 

NNODEI = Number of nodes at inlet station 

CP = Specific heat 


Total temperature 2t inlet 


tt 


ie 


a2 





Be - 


KK = Iteration counter 


NTE = Element type vector 


ALP = Nodal absolute flow angle vector 

TWEL = edad, Soe tangential velocity vector 
Ba = Nodal relative flow angle vector 

HS = Nodal static enthalpy vector 


As shown in Fig 10, the basic calculation procedure 
begins with calculating the required energy and momentun 
values at the inlet station. At this point, beginning with 
element one, the element type is interrogated to distinguish 
between duct, rotor, and stator elements. If the element is 
in a duct region, then the streamline intersections for 
local nodes 2,6,7,8 and 1 (Fig 7) are determined along with 
the associated values of energy and angular momentum. For 
the rotor and stator elements, one must initially find the 
energy and momentum values at local nodes 3,4,5 (Fig 7) due 
to the discontinuities imposed by the blade edges. Once 
these calculations are performed, then the process for the 
remaining nodes in the element proceeds in a Similar fashion 


to the duct eiements. 


After all the elements have been cycled trrough, the 
new distributions of nodal angular momentum and energy are 
returned to the main program for further computations. 
Specifically, these values will be used by the next 
subroutine , FCAL, for caiculation of the right hand _ side 


Vector. 
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Figure 10 - SUBROUTINE SLINE 
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FCAL calculates the right hand side vector as 
defined by  Gquations “{11.A.31) and (1128.21). Using tie 
identical coordinate transformations Lou nunerical 
integration as described in Section II.Cc, the final equation 
to be coded is the following, 

1 | 
f : Na 2NNe: (Lr) 
PENN Fue, LEE 2 3 


as ea (III.c. 4. 1) 
tT (2 2) Fl -H; dd Td4 dy 


where isentropic flow is assumed, and, 


W- 


; anguiar momentum 


(Til sc.0 2) 


ry 
iN 


total enthalpy 


The arguement list is defined below. FR addatton, 
only those variables in the list which have not been defined 


previously are described. 


SUBROULING FCAL(P pW ,i7 O87 CAG UVEL, RC, ZC, an, LV EG NPS 
,NODE,NN,NE,NNFSP, TWEL, NTE) 


F = Right hand side vector, f(r,Z) 

W = vector of gaussian quadrature coefficients 

ZA = Vector of 4: gaussian quadrature points 

EA = Vector of Ne gaussian quadrature points 

NFS = Vector containing nodes where the right hand sile 


a3 








is to be specified 


NN = number of nodes 

NE = number of elements 

NNPSP = Number of nodes where the right hand side is 
specified 


Fg 11 depicts the basic flowchart for the 
subroutine. To initialize the procedure, one begins with 
the first node (upper right hand corner) of the first 
element. A Switch is then applied which determines if the 
right hand side is to be calculated at the node or if a 
stream function value has been specified. This information 
1s transferred from the nain program through the arguement 
list. Once the node is allowed through the switch, then the 
integration process is started at the first integration 
POLit.. AS In? Section) Lilwa.7, “the flowchart” depices <a 
three-point Gauss Quadrature scheme. After the integration 
has been completed, a switch determines if all the local 
nodes in the element have been cycled through and if so, 
then the assembly of the elemental vector, FS, 1s performed 
to build the system right hand sid= vector , F. Finally, 
the subroutine determines if all the elements have been 
examined in order to signal completion of the right hand 
side vector. At this point, the vector, F({r,z), is returned 


to the main program for problem solution. 








BEGIN WITH lst 
EUEEN Tesi = 4 


II = II+l 
BEGIN WITH lst 
NODE OF ELEMENT. I = 1 


IS F(R,Z) TO ——_ T= 1 
SPECIFIED AT 


NEXT NODE 
THIS NODE 2 No = 











BEGIN WITH FIRST 
GAUSSIAN INTEGRATION 
POINT, 4 5:N J = 1 











FIND SHAPE FUNCTIONS, 
JACOBIAN, INVERSE OF 
JACOBIAN DET J 














PENDS NV og ly Vs 
se i’? ei 




















ZN.R, 
3 NEXT GAUSS 
FIND = (/ DW. POINT 
i=l 1 
q 
FIND # (L 4)4, NO 
i=] z 
FIND ELEMENTAL 
[NEXT NODE | RHS VECTOR  FS(TI) 
NO YES 
ASSEMBLE ae <> 
RHS VECTOR FOR 
SYSTEM. F NO 
<a> NEXT ELEMENT 
YES 
RETURN 


Figure 11 - SUBROUTINE FCAL 
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ys Subroutine vel 


[a cet A Se — == ow 


This subroutine calculates axial and radial 
velocities and also densities at each of the nodes from a 
Known strean function distribution. As noted previously in 
Section Lise. 25> Doth wWelocity “and “density “prGctiles are 
updated after obtaining the latest value of nodal strean 


fFURCELON. 


The velocity calculation proceeds from the strean 


function equations, 


3¢ 


a 


Va 


frb or (rt c.5. 1) 
| 3 
Vez 7 orb Se. (TIT .C2552) 


where 'b' is the tangential blockage factor. Since rr, ,and 


are of the following £o6ra, 


8 
i 2 % Ni 
=i 
S 
te 2 fi Ni (3 TTC. 5.5) 
{> 


then the equation for the axial velocity, Vz becomes, 


i 9 
Va = ; z aN y, (TIT.c.5.4) 
oe pn: 2 Nog Leer er 





Again, since the shape function, Le is not an 


6 1 








Implicit Lunction of "r* and %z", one MnSt use eqtation 
(II.C.1.5) to obtain the proper derivatives for computation 
of equation (TII.c.5.4). For example, CeO equation 
(TTTse55), 


aN 


3 . 
3. 2 timo i). 38 a Tia) | (Ties 5s 3) 


At this point, with equation Cit C65) 
Substituted into equation (III.C.5.4), one has the complete 
expression for the axial velocity as functions of & ur One 
proceeds similarly for expressing the radial velocity, V , 
in terms of § ane 


In order to calculate the nodal density, one uses 


the following density relation for flows in the stator and 


duct regions. 


diy (\ | ya) Fr 


(III.c.5.6) 


where Pt is the stagnation density. 


Since, 


\V's (Ve + Ve) (1+ tie’) (ite <7) 


then, 
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Since the density appears on both sides of the 
equation, the new nodal density is obtained iteratively at 


the node. 


For the relative flows in the rotor, the following 


relation for static density is used [{Ref.14]. 
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eee 
7 am tea _ Can ee wre oer (Pit 25.9) 
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Again, the solution of the nodal density is obtained 


in an iterative fashion. 


In the following arguement IASt, only those 
variables not defined in the previous subroutine 


descriptions are noted. 
SUBROUTINE VEL(NE,NN, RC, NODE,G,RG,TT,RHOT,RHON,ZC, 
Pol, nuUy, BeULNGET , UVELAVVELZRHOS TA, NEE, ALP) 
G = Ratio of specific heats 
RG = Gas Constant 
RHOT = Total density at the inlet 


RHON = work vector which contains the new nodal density 


distribution 
RHO = Nodal static density vector 
B = nodal blockage factor vector 
RHOSTA = Static density at the inlet station 
Tne “basic flowchart Zor SUBROUTING VEL. 15 Shown on 


Fig 12. Beginning with the first node of the first element, 


the Jacobian matrix (squation ({II.c.1.4)) and its inverse 
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are found. At this point the partial derivatives with 
respect to ‘'r' and 'z' of the shape functions are found as 
noted in equation (ITI.c.5.5). A Switch then allows those 
nodes not at the inlet station to pass and calculates the 
new density and velocities at the nodes. For those nodes at 
the inlet, the velocities and static densities are retained 
at the given inlet conditions. This is done to maintain 
boundary Cond.ti10n) 2nteqrity for the solution. After 
cycling through all elements, the subroutine returns the new 
nodal velocity and density distributions to the main 


program. 








BEGIN WITH lst ELEMENT, LI=1 
BEGIN WITH lst NODE OF ELEMENT 
I=1 TI=I1+1 


FIND JACOBIAN t= iol 


FIND INVERSE OF JACOBLAN 


SN: \N: 
ee eee: We 
FIND 7 1 







NO 


FIND Ve, Vz 
<r NO NEXT NODE 
YES 
ren > ae NEXT ELEMENT 
YES 


Figure 12 - SUBROUTINE VEL 
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6. Subroutine mplot 


= ee cam co se ee ce 


This subroutine utilizes the Calcomp plotter to 


depict the mesh topology of the machine under observation. 
SUBROUTINE MPLOT(RC,ZC,NODESNN NE) 
This completes the description of the main oprogran 


and associated subroutines. In the next section, a test 


case is carried through from data input to final results. 








The program was tested by using published performance 
data (Ref.12] of the NASA Task-1 stage transonic compressor. 
The compressor was discretized into twenty-eight elements 
and 107 nodes with 15 axial calculation stations (Fig 6). 
The speed was 0.5 design speed with a mass flow of 107.6 
lbm/sec. In addition, uniform flow was assumed both at the 
inlet and outlet stations. Turning angles for the rotor and 
stator were pre calculated assuming uniform conditions at 
the rotor inlet and using NASA SP-36 blade correlation data 
(Ref.13}. These absolute and relative flow angles were 
assumed constant throughout the iterative procedure as they 
were an integral part of the input data. Th2 Appendix 
contains a listing of the input data and output results for 
the NASA Task-1 transonis compressor with test conditions 
noted. To compare the accuracy of the predicted flow with 
actual laboratory observations, computed axial velocity 
profiles at the rotor inlet, rotor outlet, stator inlet, and 
stator outlet were compared with experimental results. In 


addition numerical results from Ref.7 were also compared. 


Fig 13-16 show the conputer predictions plotted with the 
experimental values and the numerical solutions obtained DY 
Hirsch and Warzee. The profiles shown were obtained after 
ten iterations and using a relaxation factor of 0.2. The 
figures show that the best overall agreement with 
experimental data occurred in the stator inlet and outlet. 
In this region the worst error was 17% which occurred at the 
Stavor..=1p inlet. The average error Chrougnouc Cie Stator 


region with respect to experimental data was 6.6%. 





The EOLOGE hub and tip outlet area exhibited 
instabilities in density convergence using equation 
Tie 25 os Specifically, the density solution converged to 
Within 8% at the rotor outlet tip and hub. It was found 
that by not allowing the nodal density at these nodes to go 
DelOW 42 Critical vatue of 0.06 lbaveuvit, the solution for 
the stream function converged. By allowing the nodal 
densities at the rotor outlet tip and hub to go below this 
critical value, the computed velocities at these nodes 
became increasingly large and the arguement within the 
brackets of equation III.C.5.9 became less than one. This 
prevented continuation of the iterations for the stream 
EPUNCtION -<SoOLUeHON. in*- addition,» he rotor » tip “outlet 
exhibited more instability than the rotor hub outlet. The 
Static density at the rotor hupd outlet oscillated about a 
value of 0.062 lom/cu ft while the rotor tip outlet was 
constantly driven to the critical value of 0.06 lbm/cu ft. 
One method attempted to alleviate this problen was the 
following. Since a half-interval iteration routine was 
used, one trial run involved reversing the direction of 
consecutive guesses when the density iteration did not 
converge. It was found however, that after three to four 
iterations of the systen of equations, the static densities 
at the rotor outlet tip and hub were again driven to smaller 
and smaller values which led to instability once more. The 
nodal densities converged at all interior points of the 
rotor edge and mid-blade regions and also at all the rotor 
inlet nodes. By including all rotor nodes, the average 


error with respect to experimental data was 27.5%. 


Fig 17 shows a plot of convergence criteria,€, versus 
the number of iterations for a relaxation factor of 0.2. 
The stability of convergence is shown to initially decrease 
and then after the third iteration oscillates about an 
approximate value of 28%. It is important to note that this 


curve represents the maximum value of € aS shown in equation 
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the curve in actuality represents 


in addition, 
function values in the 


of nodal stream 
Since in tact 


(Tie. 2. So). 
the oscillation 
rotor/stator regions 
non-linearity is the greatest. 


this is where the 





V. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDY 


Agreement with both experimental data and numerical 
solution of Ref.7 was best in the stator region to within 
8%. Predicted axial velocity profiles in the rotor inlet 
area were within 26.24 of experimental results. The 
instabilities with respect to static density solutions are 
prevalent. One of the reasons for this numerical 
disagreement with Hirsch and Warzee is the isentropic 
assumption imposed by the present program. Recommendations 
for further study on the project ainclude the addition of 
entropy variations in the rotor and stator blade regions. 
This would necessitate the use of blade correlation data 
fRef.13] for loss predictions and involve additional input 


data plus program additions to Subroutine's SLINE and FCAL. 
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Figure 14 - AXIAL PROFILE AT ROTOR OUTLET 
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Figure 16 - AXIAL PROFILE AT STATOR OUTLET 
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APPENDIX D 


CALCULATION OF ROTOR ELEMENT FLOW ANGLES 


The following 1s a brief synopsis of the procedure 
contained in Ret.13 for calculating the outlet relative 
flow angles in a rotor element from the given inlet relative 
flow angle and blade solidity. The reader is referred to 
Retr, 134 Chapter Yi, “£0 Specific details of. Bow Soeed 


COre= St 10m cata. 


AS stated 1 Section Lli~<a, wattorm £low conditions Jae 
the rotor blade edges were assumed. This assumption couvdled 
With knowledge of the mass flow rate and rotational  sveed, 
enanles one to calculates the inlet relative flow angl2,3,, 


as shown in Fig 18. 


From blade geometry information, the blade solidity,‘ , 
q = = (1) 


is obtained. At this point, P. and € are given and one may 
calcitate Pa , the rotor outlet relative flow angle from 
correlation curves depicted in Ref 13. The equation used to 


deternine Pa is the following, 


Ba= Ke + (2) 


where Ke is the angle between the tangent to the biad@ mean 





camber line and the axial direction (Fig 18). This is 
obtained from the blade geometry data. § is the low speed 
deviation angle which is obtained from the correlation 
curves in Ref. 13. The following equations show the 


relationship between § and the correlation data. 


f=. + Wo (3) 


. : (Ko). (KS), (6.),. (4) 


The variables a, K8). ,«S), and Se"), , are all values 
which are obtained from the correlation curves and are all 
functions of the given blade geometry. The quantity, $ , 1s 
the blade camber angle and again 1S obtained from the blade 
geometry data. Once all the variables are obtained fron the 
correlation data, equation (4) is solved for the deviation 
angle for an uncambered blade section, §,, and then equation 
(3) is solved for the deviation angle, S . One now 
calculates B2 from equation (2) for the blade element. With 

Ra now a known quantity, one now calculates the absolute 


flow angle, from uniform flow assumptions. 


27 

An example follows for node numbers 43 and 57 (Fig $9). 
From Ref. [12], Table iI, the following quantities are 
obtained assuming the angle of incidence, 1, (Fig 13) is 
zero and therefore the inlet relative flow angle, B, aS 


equal to Ky 


PS 61.88 
= 1.3062 
db = 6.95" 





¥, = 54.93" 


ae 
—) max = 0.035 


£18 tad10S.= 18.25 10 
hub radius = 9.125 in 
Assuming  uniforn ELOw at’ the -roteor -i1ni]et -ana a 


rotational speed of 4359.5 RPM, the following quantities are 


determined from the rotor inlet velocity diagram (Fig 18). 


: mM (107.6 lb /sec) (4a int/Fe) 
i = = 

T (0.08 \bm/et?) 1 (18.257 - 4.1257) IN? 
Vm = 246.802 Ft/sec 


where the area A, is determined from the hub and tip radii 


and the density is assumed to be 0.08 lbm/cu ft. 


Now one is ready to obtain the correlation data. From 
Ref.[13}], Fig 162, with B= 61.88° andT = 1.3062, 


5.) = D50 


From Ref.(13], Fig 162, with A = 61.388 and & =1.3062, 


f= ene oO 
From Ret... 15), Pig (17/2, with t7c) max = 0.0350, 


KM) t = 6.29 
From Ref.{13], page 222, one uses the following value of 


(K§) sh for 65-series blades, 


KSish = 129 
At this point all the necessary data has peen obtained for 


equations (3) and (4), 


C = (1.0) (0.29) (2.50) = 0.725° 


From equation (3), 


§ = 0.72574 0.235(6.95) = 2.36° 


Finally, equation (2} gives the desired value of Be 


~ 
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B = 50,95 + 2.36 2957.00" 
At this point the relative flow angle for node 57 has been 
obtained, , = 57.29°. These two values of relative flow 
angles, pi = 61.88° for node 43 and P. = 57.29" for node 57, 
acre then read in the program as input data for numerical 


computation. 


This process is repeated at each required blade element 


section for the proper outlet relative flow angle. 
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AXIAL DIRECTION, Z 





Figure 18 - NOMENCLATURE FOR CASCADE BLADE 
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